N = 10;
K = 10;

wt = zeros(N,1);
ht = zeros(N,1);

MU = [140 60];
C  = K * [20 1; 1 15];

[eigenvec, eigenval ] = eig(C);
[largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);

% Get the largest eigenvalue
largest_eigenval = max(max(eigenval));

% Get the smallest eigenvector and eigenvalue
if(largest_eigenvec_ind_c == 1)
    smallest_eigenval = max(eigenval(:,2));
    smallest_eigenvec = eigenvec(:,2);
else
    smallest_eigenval = max(eigenval(:,1));
    smallest_eigenvec = eigenvec(1,:);
end

% Calculate the angle between the x-axis and the largest eigenvector
eigenvalx=eigenval;
anglex = atan2(largest_eigenvec(2), largest_eigenvec(1));

for i=1:N
    r = mvnrnd(MU,C);
    wt(i) = round(r(1));
    ht(i) = round(r(2));
    
    fprintf ( '%d  %d  a\n', wt(i), ht(i) );
end

fprintf ( '%d  %d  b\n', round(mean(wt)), round(mean(ht)) );

MU = [60 100];
C = K * [10 -3;-3 15];

[eigenvec, eigenval ] = eig(C);
[largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);

% Get the largest eigenvalue
largest_eigenval = max(max(eigenval));

% Get the smallest eigenvector and eigenvalue
if(largest_eigenvec_ind_c == 1)
    smallest_eigenval = max(eigenval(:,2));
    smallest_eigenvec = eigenvec(:,2);
else
    smallest_eigenval = max(eigenval(:,1));
    smallest_eigenvec = eigenvec(1,:);
end

% Calculate the angle between the x-axis and the largest eigenvector
eigenvaly=eigenval;
angley = atan2(largest_eigenvec(2), largest_eigenvec(1));


for i=1:N
    r = mvnrnd(MU,C);
    wt(i) = round(r(1));
    ht(i) = round(r(2));

    fprintf ( '%d  %d  c\n', wt(i), ht(i) );
end

fprintf ( '%d  %d  d\n', round(mean(wt)), round(mean(ht)) );

MU=[160 160];
C = K * [30 20; 20 25];


[eigenvec, eigenval ] = eig(C);
[largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);

% Get the largest eigenvalue
largest_eigenval = max(max(eigenval));

% Get the smallest eigenvector and eigenvalue
if(largest_eigenvec_ind_c == 1)
    smallest_eigenval = max(eigenval(:,2));
    smallest_eigenvec = eigenvec(:,2);
else
    smallest_eigenval = max(eigenval(:,1));
    smallest_eigenvec = eigenvec(1,:);
end

eigenvalz=eigenval;
% Calculate the angle between the x-axis and the largest eigenvector
anglez = atan2(largest_eigenvec(2), largest_eigenvec(1));


for i=1:N
    r = mvnrnd(MU,C);
    wt(i) = round(r(1));
    ht(i) = round(r(2));

    fprintf ( '%d  %d  e\n', wt(i), ht(i) );
end

fprintf ( '%d  %d  f\n', round(mean(wt)), round(mean(ht)) );

sqrt(eigenvalx)
180*anglex/(2*pi)
sqrt(eigenvaly)
180*angley/(2*pi)
sqrt(eigenvalz)
180*anglez/(2*pi)